# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
## header ##
# R version 4.1.1 (2021-08-10)
# Platform: x86_64-w64-mingw32/x64 (64-bit)

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
## setup ----

rm(list=ls())

# check whether other packages than base packages are loaded
# if yes, they are detached to avoid
loadedpackages <- names(sessionInfo()$otherPkgs)
if (length(loadedpackages>0)) {
  lapply(paste("package",loadedpackages, sep=":"), detach, character.only = TRUE, unload = TRUE)
}

# define and set working directory
mywd <- "your working directory"
setwd(mywd)

# check if required folders exist
folders <- c("data","scripts","out")
all(sapply(folders, dir.exists))

# install required packages from CRAN
p_needed <- c("devtools","stargazer","survival","tidyverse","truncnorm", "zoo")
# package versions used: devtools 2.4.3; stargazer 5.2.2; survival 3.2-13; tidyverse 1.3.1; truncnorm 1.0-8; zoo 1.8-9
packages <- rownames(installed.packages())
p_to_install <- p_needed[!(p_needed %in% packages)]
if (length(p_to_install) > 0) {
  install.packages(p_to_install)
}

# load packages
lapply(p_needed, require, character.only = TRUE)

# coxed has two minor errors in its code that prevent it from working with a tvc, which we fixed in this fork
# this is a fork based on version 0.3.7
# here, I check if you already have the package installed and, if you do, will re-install the version you currently have at the end of the file
# CAUTION: if you installed from CRAN you will get the same version you had, if you installed from Github you will get the most recent version from there
# if you therefore need the old version from Github you have currently installed, do not install this
try(coxedversion <- utils::packageDescription("coxed")$Version)
install_github("imrem/coxed", upgrade="never")
require(coxed)

rm(list=setdiff(ls(), "coxedversion"))

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
## create figures 1 and 2 ----
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# this file is based on huge text files and is therefore prepared in advance
# it contains all applause per month between government parties by 10,000 words spoken by the party receiving applause
load(file="data/applause_by_10k.RDATA")

# split Merkel III away from the rest of the data so it doesn't get connected to Merkel I (same parties)
merkel3 <- applause_germany[(applause_germany$party_from_to=="CDU/CSU for SPD" & applause_germany$date>'2013-11-01')|
                              (applause_germany$party_from_to=="SPD for CDU/CSU" & applause_germany$date>'2013-11-01'),]

applause_germany <- applause_germany[!(applause_germany$party_from_to=="CDU/CSU for SPD" & applause_germany$date>'2013-11-01')&
                                       !(applause_germany$party_from_to=="SPD for CDU/CSU" & applause_germany$date>'2013-11-01'),]

fig1 <- ggplot(data=applause_germany, aes(x=date, y=applauseby10kwords, group=party_from_to)) +
  geom_smooth(aes(group=party_from_to, color=role, linetype=role), 
              span = 0.2, size=1, se=F) +
  geom_smooth(data=merkel3, aes(x=date, y=applauseby10kwords, 
                                group=party_from_to, color=role, linetype=role), span = 0.2, size=1, se=F)

fig1 <- fig1 + labs(title="", y="Applause", x="")+
  theme(axis.text.x = element_text(angle = 90, vjust = 1))+
  scale_color_manual(values=c("gray40", "black"))+
  scale_linetype_manual(values=c("dashed", "solid"))+
  ylim(0, 45) + scale_x_date(date_breaks = "1 year", date_labels = "%Y")   + 
  theme(legend.title = element_blank())

# add additional lines and text
fig1 <- fig1 + geom_vline(xintercept = as.numeric(as.Date("1998-10-01")), linetype=1) +
  geom_vline(xintercept = as.numeric(as.Date("2002-10-01")), linetype=1) +
  geom_vline(xintercept = as.numeric(as.Date("2005-11-01")), linetype=1) +
  geom_vline(xintercept = as.numeric(as.Date("2009-10-01")), linetype=1) +
  geom_vline(xintercept = as.numeric(as.Date("2013-12-17")), linetype=1) +
  theme_classic() +
  theme(axis.text = element_text(size=12), 
        axis.title = element_text(size=12),
        legend.title = element_blank()) + 
  annotate("text", x=applause_germany$date[23], y=4, label="Schr�der I", size = 5, fontface =2) +
  annotate("text", x=applause_germany$date[23], y=2, label="SPD-Greens", size = 5) +
  annotate("text", x=applause_germany$date[64], y=4, label="Schr�der II", size = 5, fontface =2) +
  annotate("text", x=applause_germany$date[64], y=2, label="SPD-Greens", size = 5) +
  annotate("text", x=applause_germany$date[102], y=4, label="Merkel I", size = 5, fontface =2) +
  annotate("text", x=applause_germany$date[102], y=2, label="CDU/CSU-SPD", size = 5) +
  annotate("text", x=applause_germany$date[145], y=4, label="Merkel II", size = 5, fontface =2) +
  annotate("text", x=applause_germany$date[145], y=2, label="CDU/CSU-FDP", size = 5) +
  annotate("text", x=merkel3$date[27], y=4, label="Merkel III", size = 5, fontface =2) +
  annotate("text", x=merkel3$date[27], y=2, label="CDU/CSU-SPD", size = 5) 

tiff("out/figure1.tif", res=600, compression = "lzw", height=7, width=14, units="in")
fig1
dev.off() 

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
### Figure 2

# split the data set into three parts - this is unfortunately necessary as Sch�ssel III and Kurz I are only in our data for very short time spans and need more smoothing to be displayed
kurz1 <- applause_austria[(applause_austria$party_from_to=="�VP for FP�" & applause_austria$date>'2010-01-01')|
                            (applause_austria$party_from_to=="FP� for �VP" & applause_austria$date>'2010-01-01'),]

schuessel3 <- applause_austria[(applause_austria$party_from_to=="�VP for BZ�" | applause_austria$party_from_to=="BZ� for �VP"),]

applause_austria <- applause_austria[!(applause_austria$party_from_to=="�VP for FP�" & applause_austria$date>'2010-01-01')&
                                       !(applause_austria$party_from_to=="FP� for �VP" & applause_austria$date>'2010-01-01')&
                                       (!(applause_austria$party_from_to=="�VP for BZ�" | applause_austria$party_from_to=="BZ� for �VP")),]

fig2 <- ggplot(data=applause_austria, aes(x=date, y=applauseby10kwords, group=party_from_to)) +
  geom_smooth(aes(group=party_from_to, color=role, linetype=role), 
              span = 0.25, size=1, se=F) +
  geom_smooth(data=kurz1, aes(x=date, y=applauseby10kwords, 
                              group=party_from_to, color=role, linetype=role), span = 0.5, size=1, se=F)+
  geom_smooth(data=schuessel3, aes(x=date, y=applauseby10kwords, 
                                   group=party_from_to, color=role, linetype=role), span = 0.7, size=1, se=F) 

fig2 <- fig2 + labs(title="", y="Applause", x="")+
  theme(axis.text.x = element_text(angle = 90, vjust = 1))+
  scale_color_manual(values=c("gray40", "black"))+
  scale_linetype_manual(values=c("dashed", "solid"))+
  ylim(0, 45) + scale_x_date(date_breaks = "1 year", date_labels = "%Y", limits = c(as.Date('2003-05-01'), as.Date('2018-08-01')))   + 
  theme(legend.title = element_blank()) 

# add additional lines and text
fig2 <- fig2 + geom_vline(xintercept = as.numeric(as.Date("2007-03-01")), linetype=1) +
  geom_vline(xintercept = as.numeric(as.Date("2008-12-02")), linetype=1) +
  geom_vline(xintercept = as.numeric(as.Date("2013-12-16")), linetype=1) +
  geom_vline(xintercept = as.numeric(as.Date("2016-05-17")), linetype=1) +
  geom_vline(xintercept = as.numeric(as.Date("2017-12-1")), linetype=1) +
  theme_classic() +
  theme(axis.text = element_text(size=12), 
        axis.title = element_text(size=12),
        legend.title = element_blank()) + 
  annotate("text", x=applause_austria$date[41], y=4, label="Sch�ssel II", size = 5, fontface =2) +
  annotate("text", x=applause_austria$date[41], y=2, label="�VP-FP� (until 2005)", size = 5) +
  annotate("text", x=applause_austria$date[41], y=0, label="�VP-BZ� (since 2005)", size = 5) +
  annotate("text", x=applause_austria$date[55.9], y=4, label="Gusenbauer", size = 5, fontface =2) +
  annotate("text", x=applause_austria$date[55.9], y=2, label="SP�-�VP", size = 5) +
  annotate("text", x=applause_austria$date[93.5], y=4, label="Faymann I", size = 5, fontface =2) +
  annotate("text", x=applause_austria$date[93.5], y=2, label="SP�-�VP", size = 5) +
  annotate("text", x=applause_austria$date[134.7], y=4, label="Faymann II", size = 5, fontface =2) +
  annotate("text", x=applause_austria$date[134.7], y=2, label="SP�-�VP", size = 5) +
  annotate("text", x=applause_austria$date[156], y=4, label="Kern", size = 5, fontface =2) +
  annotate("text", x=applause_austria$date[156], y=2, label="SP�-�VP", size = 5) +
  annotate("text", x=as.Date('2018-08-01'), y=4, label="Kurz I", size = 5, fontface =2) +
  annotate("text", x=as.Date('2018-08-01'), y=2, label="�VP-FP�", size = 5) 

tiff("out/figure2.tif", res=600, compression = "lzw", height=7, width=14, units="in")
fig2
dev.off() 

rm(list=setdiff(ls(), "coxedversion"))

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
## face validity ----
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

load("out/coalitionmood.RDATA")

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
### Figure 3

fig3 <- ggplot(mood_de) + 
  stat_smooth(aes(x = date, y = mood, group = cabinet), method = "loess", span = 0.5, size=1, se=FALSE, colour = "black") +
  stat_smooth(aes(x = date, y = conf.low, group = cabinet), method = "loess", span = 0.5, size=0, se=FALSE, colour = NA) +
  stat_smooth(aes(x = date, y = conf.high, group = cabinet), method = "loess", span = 0.5, size=0, se=FALSE, colour = NA) 

# build plot object for rendering 
gg1 <- ggplot_build(fig3)

# extract data for the loess lines from the 'data' slot
df_loess <- data.frame(x = as.Date(gg1$data[[1]]$x,origin = "1970-01-01"),
                  ymin = gg1$data[[2]]$y,
                  ymax = gg1$data[[3]]$y) 
df_loess <- df_loess %>%
  mutate(cabinet = ifelse(x<="2002-10-01", "Schr�der I",
                          ifelse(x >= "2002-10-01" & x <="2005-11-01", "Schr�der II", 
                                 ifelse(x >= "2005-12-01" & x <="2009-10-01", "Merkel I", 
                                        ifelse(x >= "2009-11-01" & x <="2013-12-01", "Merkel II", 
                                               ifelse(x >= "2014-01-01" & x <="2018-02-01", "Merkel III", 
                                                      ifelse(x >= "2018-03-01", "Merkel IV", "")))))))


# use the loess data to add the 'ribbon' to plot 
fig3 <- fig3 + geom_ribbon(data = df_loess, aes(x = x, ymin = ymin, ymax = ymax, group = cabinet),
                               fill = "grey", alpha = 0.65)

fig3 <- fig3 + labs(title="", y="", x="") +
  theme(axis.text.x = element_text(angle = 90, vjust = 1)) +
  scale_color_manual(values=c("black")) +  
  scale_x_date(date_breaks = "1 year", date_labels = "%Y") +
  scale_y_continuous(breaks=seq(0,10,1), limits = c(0, 10))

# add additional lines and text

fig3 <- fig3 + geom_vline(xintercept = as.numeric(as.Date("1998-10-01")), linetype=1) +
  geom_vline(xintercept = as.numeric(as.Date("2002-10-01")), linetype=1) +
  geom_vline(xintercept = as.numeric(as.Date("2005-11-01")), linetype=1) +
  geom_vline(xintercept = as.numeric(as.Date("2009-10-28")), linetype=1) +
  geom_vline(xintercept = as.numeric(as.Date("2013-12-01")), linetype=1) +
  geom_vline(xintercept = as.numeric(as.Date("2005-07-1")), linetype=2, color="gray60") +
  geom_vline(xintercept = as.numeric(as.Date("2012-09-13")), linetype=2, color="gray60") +
  geom_vline(xintercept = as.numeric(as.Date("2001-11-16")), linetype=2, color="gray60") +
  theme_classic() +
  theme(axis.text.y=element_text(size=12),
        axis.text = element_text(size=12), 
        axis.title = element_text(size=12),
        legend.position="none")


fig3 <- fig3 + 
  annotate("text", x=mood_de$date[22.2], y=0.7, label="Schr�der I", size = 5, fontface =2) +
  annotate("text", x=mood_de$date[22.2], y=0.35, label="SPD-B90/GR�NE", size = 5) +
  annotate("text", x=mood_de$date[64.2], y=0.7, label="Schr�der II", size = 5, fontface =2) +
  annotate("text", x=mood_de$date[64.2], y=0.35, label="SPD-B90/GR�NE", size = 5) +
  annotate("text", x=mood_de$date[101.5], y=0.7, label="Merkel I", size = 5, fontface =2) +
  annotate("text", x=mood_de$date[101.5], y=0.35, label="CDU/CSU-SPD", size = 5) +
  annotate("text", x=mood_de$date[145.5], y=0.7, label="Merkel II", size = 5, fontface =2) +
  annotate("text", x=mood_de$date[145.5], y=0.35, label="CDU/CSU-FDP", size = 5) +
  annotate("text", x=mood_de$date[191], y=0.7, label="Merkel III", size = 5, fontface =2) +
  annotate("text", x=mood_de$date[191], y=0.35, label="CDU/CSU-SPD", size = 5) +
  annotate("text", x=mood_de$date[38], y=10, label="Decision on sending", size = 4, colour="gray40", hjust = 0) +
  annotate("text", x=mood_de$date[38], y=9.7, label="troops to Afghanistan", size = 4, colour="gray40", hjust = 0) +
  annotate("text", x=mood_de$date[79], y=10, label="Vote of no", size = 4, colour="gray40", hjust = 0) +
  annotate("text", x=mood_de$date[79], y=9.7, label="confidence", size = 4, colour="gray40", hjust = 0) +
  annotate("text", x=mood_de$date[157.5], y=10, label="Final decision", size = 4, colour="gray40", hjust = 0) +
  annotate("text", x=mood_de$date[157.7], y=9.7, label="on ESM", size = 4, colour="gray40", hjust = 0)

tiff("out/figure3.tif", res=300, compression = "lzw", height=7, width=14, units="in")
fig3
dev.off() 

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
### Figure 4

fig4 <- ggplot(mood_at) + 
  stat_smooth(aes(x = date, y = mood, group = cabinet), method = "loess", span = 0.7, size=1, se=FALSE, colour = "black") +
  stat_smooth(aes(x = date, y = conf.low, group = cabinet), method = "loess", span = 0.7, size=0, se=FALSE, colour = NA) +
  stat_smooth(aes(x = date, y = conf.high, group = cabinet), method = "loess", span = 0.7, size=0, se=FALSE, colour = NA) 

# build plot object for rendering 
gg1 <- ggplot_build(fig4)

# extract data for the loess lines from the 'data' slot
df_loess <- data.frame(x = as.Date(gg1$data[[1]]$x,origin = "1970-01-01"),
                  ymin = gg1$data[[2]]$y,
                  ymax = gg1$data[[3]]$y) 
df_loess <- df_loess %>%
  mutate(cabinet = ifelse(x<="2005-05-01", "Sch�ssel II",
                          ifelse(x >= "2005-05-01" & x <="2006-12-1", "Sch�ssel III", 
                                 ifelse(x >= "2007-01-1" & x <="2008-11-1", "Gusenbauer", 
                                        ifelse(x >= "2008-12-1" & x <="2013-11-1", "Faymann I", 
                                               ifelse(x >= "2013-12-1" & x <="2016-04-1", "Faymann II", 
                                                      ifelse(x >= "2016-05-1" & x <="2017-11-1", "Kern", 
                                                             ifelse(x >= "2017-12-1", "Kurz", ""))))))))

# use the loess data to add the 'ribbon' to plot 
fig4 <- fig4 + geom_ribbon(data = df_loess, aes(x = x, ymin = ymin, ymax = ymax, group = cabinet),
                               fill = "grey", alpha = 0.65)

fig4 <- fig4 + labs(title="", y="", x="") +
  theme(axis.text.x = element_text(angle = 90, vjust = 1)) +
  scale_color_manual(values=c("black")) +  
  scale_x_date(date_breaks = "1 year", date_labels = "%Y", limits = c(as.Date('2003-05-01'), as.Date('2018-07-01'))) +
  scale_y_continuous(breaks=seq(0,10,1), limits = c(0, 10))

# add additional lines and text

fig4 <- fig4 + geom_vline(xintercept = as.numeric(as.Date("2007-01-01")), linetype=1) +
  geom_vline(xintercept = as.numeric(as.Date("2008-12-02")), linetype=1) +
  geom_vline(xintercept = as.numeric(as.Date("2013-12-01")), linetype=1) +
  geom_vline(xintercept = as.numeric(as.Date("2016-05-01")), linetype=1) +
  geom_vline(xintercept = as.numeric(as.Date("2017-12-1")), linetype=1) +
  geom_vline(xintercept = as.numeric(as.Date("2005-03-1")), linetype=2, color="gray60") +
  geom_vline(xintercept = as.numeric(as.Date("2008-07-8")), linetype=2, color="gray60") +
  geom_vline(xintercept = as.numeric(as.Date("2017-05-14")), linetype=2, color="gray60") +
  geom_vline(xintercept = as.numeric(as.Date("2015-09-15")), linetype=2, color="gray60") +
  theme_classic() +
  theme(axis.text.y=element_text(size=12),
        axis.text = element_text(size=12), 
        axis.title = element_text(size=12),
        legend.position="none")

fig4 <- fig4 + 
  annotate("text", x=mood_at$date[16.5], y=0.7, label="Sch�ssel II", size = 5, fontface =2) +
  annotate("text", x=mood_at$date[16.5], y=0.35, label="�VP-FP� (until 2005)", size = 5) +
  annotate("text", x=mood_at$date[16.5], y=0, label="�VP-BZ� (since 2005)", size = 5) +
  annotate("text", x=mood_at$date[33.8], y=0.7, label="Gusenbauer", size = 5, fontface =2) +
  annotate("text", x=mood_at$date[33.8], y=0.35, label="SP�-�VP", size = 5) +
  annotate("text", x=mood_at$date[72], y=0.7, label="Faymann I", size = 5, fontface =2) +
  annotate("text", x=mood_at$date[72], y=0.35, label="SP�-�VP", size = 5) +
  annotate("text", x=mood_at$date[112.7], y=0.7, label="Faymann II", size = 5, fontface =2) +
  annotate("text", x=mood_at$date[112.7], y=0.35, label="SP�-�VP", size = 5) +
  annotate("text", x=mood_at$date[134.8], y=0.7, label="Kern", size = 5, fontface =2) +
  annotate("text", x=mood_at$date[134.8], y=0.35, label="SP�-�VP", size = 5) +
  annotate("text", x=as.Date('2018-07-01'), y=0.7, label="Kurz I", size = 5, fontface =2) +
  annotate("text", x=as.Date('2018-07-01'), y=0.35, label="�VP-FP�", size = 5) +
  annotate("text", x=as.Date("2005-03-19"), y=10, label="FP�-BZ�-split", size = 4, colour="gray40", hjust = 0) +
  annotate("text", x=mood_at$date[40], y=10, label="Announcement of", size = 4, colour="gray40", hjust = 0) +
  annotate("text", x=mood_at$date[40], y=9.7, label="early elections", size = 4, colour="gray40", hjust = 0) +
  annotate("text", x=mood_at$date[121], y=10, label="Conflicts", size = 4, colour="gray40", hjust = 0) +
  annotate("text", x=mood_at$date[121], y=9.7, label="concerning", size = 4, colour="gray40", hjust = 0) +
  annotate("text", x=mood_at$date[121], y=9.4, label="'refugee crisis'", size = 4, colour="gray40", hjust = 0) +
  annotate("text", x=mood_at$date[139], y=10, label="Kurz takes", size = 4, colour="gray40", hjust = 0) +
  annotate("text", x=mood_at$date[139], y=9.7, label="over �VP", size = 4, colour="gray40", hjust = 0) 

tiff("out/figure4.tif", res=300, compression = "lzw", height=7, width=14, units="in")
fig4
dev.off() 

rm(list=setdiff(ls(), "coxedversion"))
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
## concurrent validity ----
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# load data 
load("out/mood.RDATA")
load("data/polls.RDATA")

# load unstandardized data which is needed for simulations
load("out/coefficients_unstand.RDATA")

coefficients_unstand <- rbind(model_ATcoefs, model_DEcoefs) %>%
  select(mood, std.error) %>%
  rename(unstand_mu = mood,
         unstand_sigma = std.error)
rm(list=c("model_ATcoefs","model_DEcoefs"))

# merge and add additional variables
### AUSTRIA
df_AT <- left_join(mood_at, df_elec_AT) %>%
  mutate(joint_diff = na.approx(joint_diff),
         spr = na.approx(spr),
         share_comp = na.approx(share_comp)) %>%
  mutate(spr_lag = lag(spr, n = 1), 
         spr_lead = lead(spr, n = 1),
         diff_lag = lag(joint_diff, n = 1),
         diff_lead = lead(joint_diff, n = 1),
         share_comp_lag = lag(share_comp, n = 1),
         share_comp_lead = lead(share_comp, n = 1),
         mood_lag = lag(mood, n = 1)) 

### Germany
df_DE <- left_join(mood_de, df_elec_DE) %>%
  mutate(joint_diff = na.approx(joint_diff),
         spr = na.approx(spr),
         share_comp = na.approx(share_comp)) %>%
  mutate(spr_lag = lag(spr, n = 1), 
         spr_lead = lead(spr, n = 1),
         diff_lag = lag(joint_diff, n = 1),
         diff_lead = lead(joint_diff, n = 1),
         share_comp_lag = lag(share_comp, n = 1),
         share_comp_lead = lead(share_comp, n = 1),
         mood_lag = lag(mood, n = 1)) 

rm(list=setdiff(ls(), c("df_AT","df_DE", "coefficients_unstand","coxedversion")))

# remove variables not needed anymore
df_AT <- df_AT[,c(1:12,22:31)]
df_DE <- df_DE[,c(1:12,22:31)]
df_pooled <- rbind(df_AT, df_DE)


# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
### Table 1

concurrent_validity <- lm(mood ~  joint_diff + share_comp + honeymoon + caretaker + 
                     mood_lag + coalition, data = df_pooled)
summary(concurrent_validity)

stargazer(concurrent_validity,
          star.char = c("*", "**", "***"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          type = "html",out="out/table1.doc")

# these are the simulations based on the estimated (bootstrapped) SEs
source("scripts/functions.R")

# add unstandardized estimates 
df_sim <- cbind(df_pooled, coefficients_unstand)

# set seed and replicate model 1000 times 
set.seed(1234)
sims_concurrent <- replicate(n = 1000, sim_concurrent(df_sim, "unstand_mu", "unstand_sigma"), simplify = FALSE)

df_concurrent <- sims_concurrent %>%
  map_df(broom::tidy) %>% 
  filter(term == "share_comp") %>%
  arrange(estimate) 

# how many model replications result in significant results 
simulations <- sims_concurrent %>%
      map_df(broom::tidy) %>%
      mutate(smaller = ifelse(p.value < 0.01, 1, 0)) %>%
      filter(term != "(Intercept)" & !startsWith(term,"coalition")) %>%
      group_by(term) %>%
      summarise(mean = mean(smaller))

simulations <- simulations[c(3,5,2,1,4),]
colnames(simulations) <- c("variables", "share_p")
simulations$share_t <- NA

simulations$share_t[1] <- sims_concurrent %>%
        map_df(broom::tidy) %>%
        filter(term == "joint_diff") %>%
        pull(statistic) %>%
        {abs(.) >  abs(as.numeric(summary(concurrent_validity)[["coefficients"]][, "t value"][2]))} %>%
        mean()

simulations$share_t[2] <- sims_concurrent %>%
        map_df(broom::tidy) %>%
        filter(term == "share_comp") %>%
        pull(statistic) %>%
        {abs(.) >  abs(as.numeric(summary(concurrent_validity)[["coefficients"]][, "t value"][3]))} %>%
        mean()

simulations$share_t[3] <- sims_concurrent %>%
        map_df(broom::tidy) %>%
        filter(term == "honeymoon") %>%
        pull(statistic) %>%
        {abs(.) >  abs(as.numeric(summary(concurrent_validity)[["coefficients"]][, "t value"][4]))} %>%
        mean()

simulations$share_t[4] <- sims_concurrent %>%
        map_df(broom::tidy) %>%
        filter(term == "caretaker") %>%
        pull(statistic) %>%
        {abs(.) >  abs(as.numeric(summary(concurrent_validity)[["coefficients"]][, "t value"][5]))} %>%
        mean()

simulations$share_t[5] <- sims_concurrent %>%
        map_df(broom::tidy) %>%
        filter(term == "mood_lag") %>%
        pull(statistic) %>%
        {abs(.) >  abs(as.numeric(summary(concurrent_validity)[["coefficients"]][, "t value"][6]))} %>%
        mean()

rm(list=setdiff(ls(), "coxedversion"))

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
## predictive validity ----
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# load unstandardized data which is needed for simulations, merge with government bill data
load("data/govbills.RDATA")
load("out/coefficients_unstand.RDATA")

model_ATcoefs$date <- substr(model_ATcoefs$date,1,7)
govbills_at_unstand <- merge(govbills_at, model_ATcoefs)

model_DEcoefs$date <- substr(model_DEcoefs$date,1,7)
govbills_de_unstand <- merge(govbills_de, model_DEcoefs)

govbills_unstand <- rbind(govbills_at_unstand,govbills_de_unstand)

coefficients_unstand <- govbills_unstand %>% select(mood_lag, std.error_lag) %>%
  rename(unstand_mu_lag = mood_lag,
         unstand_sigma_lag = std.error_lag)

rm(list=c("model_ATcoefs","govbills_at_unstand","model_DEcoefs","govbills_de_unstand","govbills_unstand"))

# load data
load("data/govbills.RDATA")
load("out/mood.RDATA")

# merge mood data with government bill data
mood_at$date <- substr(mood_at$date,1,7)
govbills_at <- merge(govbills_at, mood_at)

mood_de$date <- substr(mood_de$date,1,7)
govbills_de <- merge(govbills_de, mood_de)

govbills <- rbind(govbills_at,govbills_de)


# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
### Table 2

predictive_validity <- coxph(formula=Surv(tstart, time, status) ~ mood_lag + d_gov + d_opp +
                    monthspassed + monthspassedsq + 
                    agriculture + defence + economy + education + environment +
                    finance + foreign + health + industry + interior + justice +     
                    labour + social  + coalition, 
                  id = id,
                  data=govbills)

summary(predictive_validity)

stargazer(predictive_validity,
          star.char = c("*", "**", "***"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          type = "html",out="out/table2.doc")

# estimate duration for different levels of coalition mood
# this takes a few minutes
fig5data <- NULL
for (x in seq(0,10,1)) {
  coxed_help <- coxed(predictive_validity, method = "gam",
                      newdata = mutate(govbills, mood_lag = x),
                      id="id",
                      bootstrap = T,
                      B = 200)
  y <- cbind(coxed_help$mean,x)
  fig5data <- rbind(fig5data,y)
}

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
### Figure 5

# create a figure based on these estimates
fig5data <- data.frame(fig5data)

fig5 <- ggplot(data = govbills) +
  geom_histogram(aes(x=mood_lag, y=..count../24),bins = 21, color='gray85',fill="gray85", show.legend = FALSE, size=1) +
  geom_line(data = fig5data, aes(x=x,y=mean), color= 'black', size=1.1) +
  geom_ribbon(data = fig5data,aes(x=x,y=mean,ymin=lb, ymax=ub), linetype=2, alpha=0.4) +
  labs(title= '', x = 'Coalition Mood', y='Number of days') + theme_classic()   + scale_x_continuous(breaks=seq(0, 10, 1))


tiff("out/figure5.tif", res=600, compression = "lzw", height=7, width=7, units="in")
fig5
dev.off() 


# estimate how much longer a bill takes to be accepted going from the first to the third quartile of coalition mood
coxed_help <- coxed(predictive_validity, method = "gam",
                    newdata = mutate(govbills, mood_lag = summary(govbills$mood_lag)[2]),
                    newdata2 = mutate(govbills, mood_lag = summary(govbills$mood_lag)[5]),
                    id = id,
                    bootstrap = T,
                    B = 200)

coxed_help$mean1[1]
coxed_help$mean2[1]
coxed_help$mean1[1]-coxed_help$mean2[1]

## these are the simulations based on the estimated (bootstrapped) SEs
source("scripts/functions.R")

# add unstandardized estimates 
df_sim <- cbind(govbills, coefficients_unstand)

# check function
sim_predictive(df_sim, "unstand_mu_lag", "unstand_sigma_lag")

# set seed and replicate model 1000 times 
# this takes a few minutes
set.seed(1234)
sims_predictive <- replicate(n = 1000, sim_predictive(df_sim, "unstand_mu_lag", "unstand_sigma_lag"), simplify = FALSE)

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
### Figure 6

# extract distribution of estimate of interest and plot it
df_predictive <- sims_predictive %>%
  map_df(broom::tidy) %>% 
  filter(term == "mood_sim_lag") %>%
  arrange(estimate) %>%
  mutate(upper_95 = estimate + 1.96 * std.error,
         lower_95 = estimate - 1.96 * std.error,
         upper_99 = estimate + 2.58 * std.error,
         lower_99 = estimate - 2.58 * std.error,
         sim = row_number()) 

fig6 <- ggplot(data = df_predictive, aes(y = sim, x = estimate)) +
  geom_linerange(aes(y = sim, xmin = lower_99, xmax = upper_99), color = "lightgrey", size = 0.5, alpha = .5) +
  geom_point(size = 0.75) +
  geom_vline(xintercept = 0.107, linetype = "dashed") +
  ylab("Simulation") + xlab("Coalition mood (t-1)") +
  theme_classic()

tiff("out/figure6.tif", res=600, compression = "lzw",  height=4.25, width=14, units="cm")
fig6
dev.off()

sims_predictive %>% 
  map_df(broom::tidy) %>% 
  filter(term == "mood_sim_lag") %>%
  ggplot(aes(x = statistic) ) +
  geom_density(fill = "grey", alpha = .5) +
  geom_vline(xintercept = 0)

sims_predictive %>% 
  map_df(broom::tidy) %>% 
  filter(term == "mood_sim_lag") %>%
  ggplot(aes(x = estimate) ) +
  geom_density(fill = "grey", alpha = .5) +
  geom_vline(xintercept = 0)

sims_predictive %>%
  map_df(broom::tidy) %>%
  filter(term == "mood_sim_lag") %>%
  pull(p.value) %>%
  {. <  0.01} %>%
  mean()

# the forked version of coxed is uninstalled here
detach("package:coxed")
remove.packages("coxed")

# if the package was installed before this file was executed, it is reinstalled
# if it was originally downloaded from CRAN the same version as before is installed
# if it was originally downloaded from Github the current Github version is installed
if (exists("coxedversion")==TRUE){
  tryinstall <- try(install_version("coxed", version = coxedversion, 
                                    repos = "http://cran.us.r-project.org", upgrade="never"))
  
  if (inherits(tryinstall, "try-error")==TRUE) {
    install_github("jkropko/coxed", upgrade="never")
  }
  
}
